Vibrato-ify It!!

Chris Tralie

import numpy as np
import matplotlib.pyplot as plt
import IPython.display as ipd
import librosa

Sinusoid Vibrato

Repeat at frequency $f$ that cycles $r$ times per second, at a deviation of $v$ hz up and down

$\frac{df}{dt} f(t) = f + v \cos(2 \pi r t)$

$f(t) = ft + \frac{v}{2 \pi r} \sin(2 \pi r t)$

sr = 44100
f = 440*2**(-9/12)
v = 10
r = 5
t = np.arange(sr*4)/sr # In units of seconds(!)
y = np.cos(2*np.pi*(f*t + (v/(2*np.pi*r))*np.sin(2*np.pi*r*t)))
ipd.Audio(y, rate=sr)

General Parameterized Vibrato

y, sr = librosa.load("f3_long_straight_a.wav", sr=sr)
ipd.Audio(y, rate=sr)

Given $y(t)$

Re-parameterize as $y(u(t))$

Suppose as an example that $u(t) \propto t^2$ so that it keeps playing faster and faster over time

t = np.arange(y.size)/sr
tmax = np.max(t)
ut = tmax*(t/tmax)**2
plt.plot(t, ut)
plt.xlabel("t")
plt.ylabel("u(t)")

Text(0, 0.5, 'u(t)')

ywarp = np.interp(ut, t, y)
ipd.Audio(ywarp, rate=sr)

Now let's compose the audio with a function that slows down and speeds up in a periodic fashion. This gets us a vibrato sound!

f = 440*2**(-9/12)
b = 5
v = 10
a = v/(2*np.pi*b*f)
ut = t + a*np.cos(2*np.pi*b*t)
ywarp = np.interp(ut, t, y)
ipd.Audio(ywarp, rate=sr)

The math below shows what this function composition of a sinusoid $y(t) = A \cos(2 \pi f t + \phi)$ with the above function. We see it behaves the same way as the original sinusoidal vibration equation! The only thing is that the amount of vibrato is directly proportional to the frequency. So if we have a sound with multiple harmonics, the higher harmonics will cycle more

$y(t) = A \cos(2 \pi f t + \phi)$

$y(t + a \cos(2 \pi b t)) = A \cos(2 \pi f (t + a \cos(2 \pi b t)) + \phi) $

Original: $\frac{df}{dt} f(t) = f + v \cos(2 \pi r t)$

New: $ \frac{d f(t)}{dt} = f - 2 \pi a b f \sin(2 \pi b t) $

$a = v / (2 \pi b f)$

y, sr = librosa.load("westminster.mp3", sr=sr)
ipd.Audio(y, rate=sr)
t = np.arange(y.size)/sr

f = 440*2**(-9/12)
b = 10
v = 10
a = v/(2*np.pi*b*f)
ut = t + a*np.cos(2*np.pi*b*t)
ywarp = np.interp(ut, t, y)
ipd.Audio(ywarp, rate=sr)